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Abstract 

First-principles density- functional theory has been used to investigate equilibrium geometries, 
total energies, and diffusion barriers for H as an interstitial impurity absorbed in a-Fe. Internal 
strains/stresses upon hydrogen absorption are a crucial factor to understand preferred absorption 
sites and diffusion. For high concentrations, H absorbs near the octahedral site favoring a large 
tetragonal distortion of the BCC lattice. For low concentration, H absorbs near the tetrahedral site 
minimizing the clastic energy stored on nearby cells. Diffusion paths depend on the concentration 
regime too; hydrogen diffuses about ten times faster in the distorted BCT lattice. External stresses 
of several GPa modify barriers by 10%, and diffusion rates by 30%. 

PACS numbers: 66.30.J-,68.43.Bc,71.15.Mb 



INTRODUCTION AND METHODOLOGY 



Hydrogen embrittlement is believed to be one of the main reasons for cracking of steel 
structures under stress [TH3]. High strength steels often include a ferritic core made of a- 
iron (body centered cubic lattice, BCC). To control and prevent the cracking of steel it 
is necessary to understand the chemical and physical properties of hydrogen inside BCC 
iron. In particular, it has been argued in the literature that the size of the interstitial H 
(e.g. ~ 0.3 — 0.4 A[4j) is too large to fit and diffuse easily in the iron BCC lattice[5]. This 
argument has also been used to favor absorption on the T-site because it is larger than the 
0-site[n]. Arguments based in 'apparent' sizes of atoms and a static lattice yield an useful 
basic approximation but leave out important factors. Only a truly microscopic understand- 
ing of hydrogen dissolved in iron, including diffusion processes and fully accounting for forces 
and stresses inside the host matrix could lead to answers to the many questions arising from 
experiments. Internal stresses are therefore key to understand the preferred absorption site. 
Furthermore, external stresses can be used to modify diffusion barriers because can help 
to stabilize/de-stabilize strains originated inside the material to accommodate the absorbed 
impurity. Owing to its interest, hydrogen absorbed in iron has been studied from different 
approaches, both from a theoretical and experimental point of view. Unfortunately, a clear 
consensus about fundamental questions, like the nature of the equilibrium absorption site, 
the reasons for H to prefer some regions over others, or how to modify the diffusion barriers 
has not been reached yet. Many theoretical papers have favored hydrogen absorbed in the 
tetrahedral site (T-site) [7I4T2] . some have preferred the octahedral one (0-site)[T3], others 
have reported that they are almost equivalent [S] (Fig. [T]). From an experimental point of 
view, it has been generally believed that hydrogen has low solubility and high mobility in 
bcc-iron. However, some recent experiments have attained high concentrations [151 IIS]) cit 
least locally. Furthermore, experimental evidence has been reported showing that hydrogen 
might at least partly be found on the 0-site[T2]. Finally, diffusion barriers that change with 
the amount of H admitted into the material have been reported in the literature^. In 
our calculations, we find that both sites are indeed quite similar in energies; the delicate 
energetic balance between them depending on the stress distribution created by the inter- 
stitials. In turn, this depends on the hydrogen concentration, that can be related to locally 
inhomogeneous regions or more generally to its homogeneous averaged value. 
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FIG. 1: High symmetry sites for H absorption inside a-iron: the octahedral (0-site) and tetrahedral 
(T-site) sites (grey) are shown along with iron atoms in the BCC lattice (black). 

Our methodological approach simulates the atomic absorption and diffusion of hydro- 
gen in the body centered cubic iron lattice from first-principles. Absorption in high- 
symmetry sites for different hydrogen loads (densities) has been explored and the internal 
strains/stresses have been calculated. Diffusion barriers between these stationary absorp- 
tion sites have been obtained, and the effect of external stresses on diffusion coefficients has 
been analyzed. Our ab-initio calculations are based on the Density Functional Formalism 
(DFT)[TH] and the theory of pseudopotentialsP^ . A periodic supercell is considered to 
minimize absorbate-absorbate interactions (periodicity is n x n x n with respect to a single 
cubic BCC unit cell -UC-). One H atom has been located in this supercell, determining the 
ratio H:Fe from ^ for the 1 x 1 x 1 to ^ for the 3x3x3. Within this theoretical frame- 
work the accuracy of calculations mainly depend on the quality of the pseudopotential, the 
energy cutoff (-Ect/), and the density of the k-points mesh in the irreducible part of the Bril- 
louin zone (a Monkhorst-Pack set of special k-points [20]). Calculations are performed with 
CASTEP[2I1I22], where a plane- wave basis is used to expand electronic states (spin-polarized 
bands are considered to account for magnetism in iron). An ultrasoft pseudo-potential (3d^ 
4s^) including a nonlinear core correction is used to describe iron atoms. The exchange and 
correlation contribution to the total energy is computed within a generalized gradients ap- 
proximation given by Perdew, Burke and Ernzerhof |23] . An all-bands variational method is 
used to achieve self-consistency [24j ( to improve convergence a smearing width of 0.1 eV has 



been used). Minimum thresholds to consider the calculations converged after optimization 
of the geometry are: variations in the total energy < 10~^ eV, maximum residual force on 
any atom < 0.01 eV, maximum residual stress on the UC < 0.02 GPa, maximum change 
in any atom position < 0.001 A, and maximum change in any cell vector < 0.01 A. These 
thresholds are enough for our purposes, but on some particular cases lower thresholds were 
considered when higher precision was desirable. 

REFERENCE CALCULATIONS FOR a-IRON AND H 

The adequacy of the theoretical framework and the different convergence thresholds has 
been tested by computing a number of known properties for clean a-iron and atomic and 
molecular hydrogen separately. For iron, a global energy minimum is found for a BCC lattice 
(hereafter referred as the 1x1x1) with a = b = c = 2.815 A (a = /3 = 7 = 90°), a predicted 
magnetization of 2.25 fiB, a bulk modulus of 175 GPa, and a Poisson ratio of 0.30. Residual 
forces on the atoms in the BCC unit cell are < 4 x 10~® eV/A, residual stress is < 9 x 10~^ 
GPa, and the total energy has been converged to < 10^^ eV). The computed magnitudes 
are to be compared with the experimental values 2.867 A, 2.22 jjB, 170 GPa, and 0.29 125]. 
Fractional errors are 1.8%, 1.4%, 3%, and 4% respectively, showing the accuracy of this 
kind of theoretical modeling to describe relevant properties of BCC iron. It is interesting, 
however, to notice that these quantities show oscillatory convergence as a function of the 
energy cutoff (£^ct/) and the number of k-points; convergence of physical properties has to be 
established carefully to make sure the oscillations are damped enough as to avoid agreements 
that might not be stable (Fig. [2]). Properly converged values were obtained with Edf = 375 
eV and a 25 x 25 x 25 Monkhorst-Pack (MP) mesh (455 points in the irreducible part of the 
Brillouin zone). However, values derived with Edf = 375 eV and a 14 x 14 x 14 mesh can 
be considered accurate enough for our purposes and most of our calculations below for the 
1x1x1 BCC unit cell have been performed with these parameters. Calculations in nxnxn 
unit cells (n=2, 3) have been performed with 4x4x4 and 1x1x1 MP meshes respectively. 
These theoretical calculations yield a detailed microscopic picture for the chemistry and the 
physics governing the material: iron atoms form metallic bonds of length 2.44 A, favoring a 
ferromagnetic alignment of spins mostly residing in the regions around the atoms. Bonding 
is such that the atomic s level of iron is depopulated in favor of d and p-levels (a MuUiken 
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FIG. 2: Convergence of primitive unit cell parameter ( d in A) -upper panel-, magnetization (in 
units of the Bohr magneton) -middle panel-, and bulk modulus (GPa) -lower panel- versus the 
energy cut-off (circles: 330 eV, squares: 350 eV, triangles: 370 eV), and the K x K x K index for 
the Monkhorst-Pack mesh {K = 25 results in a mesh of 455 points sampling the Brillouin zone). 

analysis shows that ~ 1.3 e is transferred). Although this work is not concerned with 
face-centre cubic iron (7- iron), we notice that this formalism predicts that the next stable 
phase of iron (FCC) is about 0.1 eV/atom higher in energy, with a = b = c = 3.479 A 
(experimental value is 3.647 A, a 5% discrepancy) , and a magnetization lower by about a 
factor of 2 due to disorder in the spin orientations. Calculations for H and H2 have been 
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performed on a 10 x 10 x 10 A super-cell (only the T point in the BZ), predicting a bond 
length and a binding energy for the molecule within 0.5% and 5% from experimental values. 
The interaction between isolated atoms of H and Fe is moderately weak, with a bonding 
distance of 1.54 A and a binding energy of Eb = Ep^u — {Epe + Eh) = —2AleV (notice 
that the binding energy is quoted from the electronic energy alone and zero point vibrations 
have not been considered). This is corroborated by the very weak interaction of H with 
Fe2; the binding energy stays nearly the same (—2.43 eV) in a geometrical configuration 
where H prefers to interact with only one of the Fe in the dimer. Molecular H2 has been 
reported to dissociate on iron surfaces [12], in addition we do not find tendency for hydrogen 
atoms inside iron to re-form molecules. As an example, we compare for the 2x2x2 unit 
cell (Fei6H2) the total energies of two H atoms located in nearby two different octahedral 
sites (H-H distance is ~ 2 A) and a single H2 molecule. The dissociated state is favored 
by E2H — Eh2 = —3.3 eV (we notice that the energy gain comes partially from the smaller 
stress that the two separated H atoms introduce in the iron lattice compared with the one 
originated from the single H2 molecule). Therefore, we conclude that it is not necessary to 
study processes related to molecular H2 inside a-iron. 

HYDROGEN IN IRON: HIGH CONCENTRATION 

Even if the density of hydrogen absorbed in iron is expected to be small on the aver- 
age, it can accumulate locally on some regions where the density can be high. Indeed, 
important changes in the elastic properties of iron are found in these regions. Furthermore, 
experimental evidence of high density phases have been reported in the literature [151 116] . 
Therefore, we find it interesting to analyze the stoichiometry Fe2H. In the BCC lattice two 
high-symmetry sites (0-site and T-site) compete to locate the interstitial atom, and diffusion 
through the iron matrix takes place along paths joining them. We simultaneously optimize 
all the geometrical parameters to reach a global minimum in the total energy and small 
residual forces in any atom in the model. The three vectors defining the periodic lattice 
have been kept orthogonal (a = /3 = 7 = 90°). The lengths of these vectors have been 
optimized taking into account the following cases: (a) independent optimization (resulting 
in tetragonal systems), (b) optimization subject to the condition a = b = c (cubic system), 
and (c) no optimization, merely scaling the clean BCC unit cell (the system is under stress. 
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specially for high concentrations of the interstitial). Starting with model (a) and allowing 
the lattice to relax completely, we have reached a best configuration with residual forces 
< 10~^ eV/A and residual stresses < 5 x 10^"^ GPa (the energy change was < 2 x 10~^ 
eV and the maximum final displacement of any atom was < 5 x 10^^ A). Under these 
conditions, the 0-site is the global minimum, the T-site being higher in energy by 53 meV 
(Table [n]). To reach this configuration, the lattice suffers a large body centered tetragonal 
(BCT) distortion, keeping the P4/MMM symmetry group (16 elements), while the optimum 
geometry for the fully relaxed lattice for H in the T-site increases the symmetry from C2 
(4 elements) to P-4M2 (8 elements). Our preferred absorption site is different from the one 
recently obtained by Jiang and Carter for the same density using a similar formalism [12j. 
There are several possible sources for this discrepancy, e.g. the pseudopotentials, the ex- 
change and correlation functional, or even the fine details of how different codes reach a 
self-consistent solution and compute forces on the atoms. We believe these are all unlikely 
sources to explain the different final configuration due to the excellent performance when 
comparing different physical magnitudes with experimental data or just between different 
theoretical methods of similar quality. Jiang and Carter predict for Fe2H that the T-site 
is preferred over the 0-one by 10 meV, while we have found the reversed conclusion by 53 
meV. The overall 63 meV discrepancy is indeed a small one, but it seems to be above the 
noise level expected from the convergence tests. The reason for the difference is most likely 
related to the final optimum configuration obtained in each calculation (our final config- 
uration is given in Table |l]). The energy landscape is a complex one, and the algorithm 
searching for a global minimum could be trapped in nearby local ones. We have made an 
effort to minimize residual forces and stresses as much as possible. To search for a global 
minimum we have used the Broyden-Fletcher-Goldfarb-Shanno algorithm that is known to 
be a robust one and has been successfully applied both in electronic structure calculations 
and in Low-Energy Electron Diffraction optimization of geometries |2SH2S]- The resulting 
BCT distortion for H in the 0-site is indeed quite big and might be hidden in a difficult 
landscape, we notice that our optimization under the restriction of an isotropic deformation 
of the unit cell {a = b = c) results in the T-site as the preferred one by 452 meV. Finally, 
we should comment that although the BCT distortion seems significant, in particular be- 
cause it would favor the BCC to FCC transformation, the discussion here is based solely 
on electronic energies, and other factors might play an important role too. Configurational 
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disorder is expected to favor T-sites. The relative contribution of entropy to the Helmholtz 
free energy for T-site and 0-site can be estimated from the simple consideration that there 
are four times more configurations for T-sites. Therefore, T-sites are favored on the grounds 
of entropic arguments by ksT log A, approximately 34 meV at room temperature, which 
tends to make the energetic more similar but still is not enough to make the T-site the 
global minimum. Furthermore, zero-point energy corrections (ZPE) may affect differently 
the energy of two absorption sites with different coordination and geometrical parameters. 
A careful determination of the phonon spectra for this system demands the use of a con- 
siderable amount of computing time and is planned for a future paper. Our own initial 
estimates, based in sketchy models that can be solved quickly, results in the 0-sitc getting 
further stabilized with respect to the T-site by ~ 20 meV. A word of caution is in order 
here since these values may have a large error bar that we have not properly estimated. 
The model used to deduce these values is as follows: the cut-off energy has been reduced to 
300 eV, the k-points mesh to a mere 6x6x5, phonons have been estimated only in the F 
point using a finite displacement method (0.00529 A) and a 2 x 2 x 2 super-cell. Therefore, 
these corrections are only to be taken as rough order of magnitude estimates, but in our 
calculations ZPE does not switch the preferred absorption site. Frequencies in F for the O 
and T sites are all real, marking the character of global and local minima for these sites. 

From a physical point of view, in this high concentration regime, the main characteristic 
of the tetragonal distortion necessary to accommodate H in the 0-site is that it is both 
large and anisotropic. First, we notice that only occupation of similarly oriented 0-sites 
over many unit cells would allow this kind of distortion. Second, fractional changes in the 
UC volume are large, ^ = 0.12, and 0.16, for 0-site and T-site respectively; this is the trace 
of the strain tensor and measures the global stress produced by the interstitial. Therefore, 
it is necessary to understand the stress distribution created around the distorted cell. Total 
energies for a model where the UC is allowed to relax keeping the cubic UC {a — b — c) 
show that when H is in O the system gains a fair amount of energy from the anisotropic 
deformation (510 meV), while when H is near T-site the gain is smaller (9 meV). Under this 
cubic restriction, the fractional change of the UC volume are 0.19 and 0.16 for 0-site and 
T-site respectively; the interstitial fits poorly in O in absence of the tetragonal distortion 
(internal pressure upon injection of the interstitial H in the clean iron BCC lattice amounts 
to several tens of GPa). From a chemical point of view, the role of the interstitial atom 
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is to debilitate the Fe-Fe nearest-neighbors (NN) interaction; Fe-Fe bond length increases 
by 0.15 and 0.08 A for H in 0-site and T-site respectively. In both sites, hydrogen tends 
to get negatively charged (0.3e), while its spin is delocalized contributing little to the total 
integrated spin density that only increases by ^ 0.15 fiB with respecto to clean BCC Fe. 
It is worth to notice that the tetragonal distortion of the UC for H absorbed on the 0-site 

= 1.38) is not far away from the distortion predicted for Bain's mechanism to transform 
from a-iron to 7-iron = 1.41) [29]. 

To compute diffusion barriers we have used a linear synchronous-transit method followed 
by quadratic refinement j30] . In addition, when a reduced number of parameters could be 
established as the more relevant ones exhaustive searches in regions of interest have been 
used to display in more detail the topology around the reaction pathway (Fig. [3]). We have 
computed diffusion barriers between the optimum configurations reached for H at 0-site and 
T-site. For the high density case, the transition state is located in the line joining 0-site 
and T-site, about half-way, with a barrier of -B = 70 meV (Fig. |4]). Thermal diffusion can be 
estimated by considering a typical vibrational frequency for H in the lattice: W ~ 10^^ Hz. 
Diffusion is quick, hopping from 0-site to T-site happens Dt{B) = 10^^ times per second 
at room temperature. Quantum diffusion via tunneling through the barrier is important for 
such a light atom as H. The tunneling rate can be estimated from a one-dimensional WKB 
approximation, DQ(B,d) = where d is the distance between minima (0.6 A). 

Dq is only lower than the thermal one, Dt, by a factor of 6. At low temperature, the 
tunneling process becomes the main mechanism, being equal in importance to the thermal 
one at about T ^ 183 K, near the liquid Nitrogen temperature. The effect of external 
stresses on the barrier is shown in Fig. |4| As discussed above, ZPE can affect the energies of 
different configurations, and in particular it can lower the already small diffusion barriers. 
The model described above, results in a lowering of the barrier as seen from the 0-site of 
^ 7 meV, and a lowering of the barrier as seen from the T-site of ~ 27 meV. We notice 
that frequencies computed for the TS now include a single imaginary value, indicating that 
this is a first order saddle point. Again, this correction should be taken as a crude estimate 
and we notice that similar calculations give a correction for the barrier of 46 meV|12]. 
We estimate that barriers corrected by ZPE should be lower than the ones obtained from 
total energies only by ~ 40 — 50%. A tensile stress of 2 GPa increases or decreases the 
barrier by about 9%; changing the thermal diffusion by about 30% at room temperature. 
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TABLE I: Optimum geometrical configurations for Fe2H near the octahedral and tetrahedral sites. 
Atomic positions, x, y, z, are given in fractional coordinates (a = /5 = 7 = 90°). Total energies (in 
eV) correspond to the optimum geometry. Fmax and Smax are the maximum residual forces on 
atoms and stresses on the cell. The charge transfer, A{q), is given in units of the electron charge, 
and the spin is given in units of hb (Mulliken analysis). 

OCTAHEDRAL 



energy (eV) spin Fmax (eV/A) Smax (GPa) 
-1748.150 2.40 3 x 10"^ 4 x lO"^ 



a (A) b (A) c (A) c/a V (A3) AV/V (%) 
2.628 2.628 3.618 1.38 24.99 0.12 



atom X y z A(g) spin 
H 0.5 0.5 0.0 -0.30 -0.02 
Fca 0.0 0.0 0.0 0.19 1.11 
Feb 0.5 0.5 0.5 0.10 1.29 

TETRAHEDRAL 



energy (eV) spin Fmax (eV/A) Smax (GPa) 
-1748.097 2.46 1 x 10"^ 5 x 10"^ 



a (A) b (A) c (A) c/a V (A3) AV/V (%) 
2.896 2.984 2.984 0.97 25.78 0.16 



atom X y z A{q) spin 
H 0.2641 0.5 0.0 -0.34 -0.02 
Fea 0.0 0.0 0.0 0.17 1.24 
Feb 0.5283 0.5 0.5 0.17 1.24 
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FIG. 3: Total energy (meV) landscapes for (a) high and (b) low concentrations cases (top and 
bottom panels), x and y are given in fractional units with respecto to the 1x1x1 BCC unit 
cell. 0-site is located at the top right corner ^,0}, and T-site at the down right and top left 
corners. '^^^ region around the Fe at the origin (higher values than 200 meV) has been 

artificially put to a constant value to improve visibility of the relevant paths from 0-site to T-site 
(a) and from T-site to T-site (b). 
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while the quantum mechanism is affected by 25%. Barriers for this high concentration 
scenario have been computed using periodic boundary conditions on a 1 x 1 timesl unit cell, 
representing diffusion of H and all the adjacent periodic copies simultaneously (notice that 
periodic boundary conditions respect the Fe2H stoichiometry at all times during the diffusion 
process). Finally, we check the internal consistency of these high-density simulations by 
considering a 2 x 2 x 2 unit cell, with the same Fe2H stoichiometry. In this unit cell there 
are 36 0-sites and 144 T-sites that can be occupied by the 8 H. This makes 10^ possibihties 
for octahedral occupation, 10^^ for tetrahedral occupation, and 10^^ for an arbitrary mixture 
of both. A proper study of this system should try to sample these configurations, which 
is beyond the scope of this work. However, a check for the consistency of our previous 
results only needs to consider a few of these configurations. To ascertain the stability 
of the energy ordering between O and T-sites we consider two configurations where the 
octahedral/tetrahedral sites are occupied accordingly with results obtained in the 1x1x1 
unit cell (including the respective tetragonal distortions) . The length of the vectors and all 
the atoms in the 2x2x2 unit cell are allowed to relax in x, y and z directions including 
72 degrees of freedom (the center of mass is constrained). A stationary point is considered 
converged when forces and stresses are below 0.002 eV/A and 0.005 GPa respectively. The 
0-site converges to a better minimum than the T-site by 0.470 eV, which is consistent with 
the value derived from the 1x1x1 simulation (0.06x8) with an accuracy of 1 .2 meV / ( 1 x 1 x 1 
unit cell). In this model, if we try to compute the barrier for a single H hopping from the 
global minimum in the 0-site to the T-site, we find that the effect of the occupation of the 
other 7 sites in octahedral positions is to transform this single T-site into a high order saddle 
point. H in this T-site relaxes spontaneously to the neighboring 0-site. Barriers computed 
with the help of the 1x1x1 periodic model correspond to transforming large regions of 
all 0-occupancy into large regions of all T-occupancy. We haven't investigated diffusion 
paths for H in the presence of a given mixture of 0-sites and T-sites due to its combinatorial 
complexity. However, this example shows the interest to further study n xn xn cells in the 
high density regime in the future, and highlights the impact of the distortion of the lattice 
on the absorption energy on single local sites. 
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FIG. 4: Diffusion barriers between T-site (j) and 0-site (j) along the actual diffusion path (see 
Fig.jsj the reaction coordinate has been scaled between both endpoints) for fully relaxed Fe2H under 
different external stresses. Triangles (a) and diamonds (d) correspond to hydrostatic pressures of 
±2 GPa (tensile and compressive respectively), squares (c) to a tensile uniaxial stress of —1 GPa, 
and circles (b) to GPa. The origin of energies has always been taken as the optimum equilibrium 
value (for this concentration, the 0-site). 

HYDROGEN IN IRON: LOW CONCENTRATION AND FINITE ELEMENTS SIM- 
ULATION 

We simulate the lower density scenario by placing one H in the 2x2x2 and 3x3x3 
unit cells (FeigH and Fe54H). At these concentrations, hydrogen is free to move around the 
lattice without interacting with other hydrogen atoms (the remaining H-H interaction on 
the 1x1x1,^ 0.06 eV, now becomes negligible). Other convergence parameters are kept 
the same, except for the k-points mesh, that has been reduced accordingly to the bigger size 
of the super-cell. Now, the global total energy minimum still happens for a BCT lattice, but 
it moves from the 0-site to the T-site, that becomes lower in energy w.r.t O by 119 and 105 
meV for FcigH and Fe54H, respectively. The change in the unit cell volume and the tetragonal 
distortion of the lattice gets smaller with decreasing density, and nearly independent of the 
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TABLE II: For different concentrations and systems, energy differences between the 0-site and 

Ay 

Vb 



T-site sites (AE in meV), and fractional changes in the unit cell volume for 0-site and T-site (^^^ 



DENSITY 

SYSTEM AE ^ lo ^ |t 



Fe2H 



tetragonal -53 0.12 0.16 
cubic 452 0.19 0.16 
1 X 1 X 1 612 0.00 0.00 



FeieH 



tetragonal 119 0.03 0.03 
cubic 169 0.03 0.03 
2 X 2 X 2 117 0.00 0.00 



tetragonal 105 0.004 0.003 
cubic 115 0.004 0.004 
3 X 3 X 3 108 0.000 0.000 



Fe54H 



absorption site, ^ = 0.03 and 0.004, and ^ = 1.04 and 1.003 respectively (Table II). The 
same calculation for a relaxed BCC lattice {a = b = c = 2.843 A) yields a difference of 169 
meV in favor of T-site. Diffusion happens now from T-site to T-site with a barrier of 82 
meV (BCT) or 127 meV (BCC) (Fig. |3](b)). This barrier implies a diffusion rate ~ 10 times 
slower than for the high density case. Tunneling is even more reduced, because the distance 
between minima is larger (^ 0.9 A). In T-site, H moderately interacts with four NN at 1.65 
A, and gets negatively charged by 0.33e. In O, H interacts mainly with two NN, at 1.57 
A, while the nearest next-neighbors (NNN) remain at 1.98 A (H gets a similar amount of 
charge, 0.32e). Typical displacements for NN w.r.t. their positions in the clean iron lattice 
are ~ 0.1 A (NN), while NNN and more distant atoms move less than ^ 0.02 A. The spin 
on the different atoms does not show changes worthwhile to comment. 

Strains to accommodate the interstitial H create a stress distribution that has been an- 
alyzed by a finite elements (FE) technique (input parameters have been obtained from our 
calculations above). Assuming that deformations are elastic, we have computed the stress 
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TABLE III: Elastic energy distribution (meV), and fractional contribution to the total one (%) 
calculated by finite elements over n x n x n regions (n = 2, 3, 4, 5). Residual energy beyond n = 3 
is below a 3% of the total one. 



2x 3x 4x 5x 

0-site 

E( meV) 195 215 219 221 
% 88 9 2 0.8 



2x 3x 4x 5x 

T-site 

E( meV) 23 24 24 24 
% 95 4 0.7 0.3 

distribution in the iron matrix and the free energy variation when H is absorbed in the 
0-site and T-site. Starting from strains obtained by DFT on the 1x1x1 super-cell, the 
FE simulation allows us to study its effect in a surrounding region of size n x n x n around 
the original volume (n = 2,5). Strictly speaking, this procedure implies introducing a de- 
formation field that has been computed using periodic boundary conditions on a bulk-like 
1x1x1 cell into an infinite medium with the same elastic constants but different boundary 
conditions (non-periodic). This is therefore a slightly inconsistent procedure that can only 
be justified 'a posteriori' by the agreement between the energies obtained from the ab-initio 
calculations and the FE approach (see below). 

All the FE simulations have been performed with the program Comsol Multiphysics|31] . 
using 0.15 A wide triangular elements. Fig. [s] gives the diagonal components of the stress 
tensor for H in 0-site and T-site (the origin of distances is located at the surface of the 
deformed 1x1x1 region). Both cases show very little residual stress beyond 6 A (which 
corresponds to the 3x3x3 supercell). The elastic energy distribution associated to these 



stresses is well converged at these distances (Table III); contributions from larger regions 
only amount to 1 or 2% of the total elastic energy. This result shows why it is not necessary 
to go beyond the 3x3x3 supercell in our ab-initio simulations. The stress distributions 
are very different for H in 0-site and in T-site: the former is very anisotropic {x ^ y ^ z) 
and its maximum value is ~ 2 times the maximum value for T-site. The stored elastic 
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energy, quadratic on the stress tensor, is approximately ten times higher for H in the 0-site 



(Table III). Therefore, storing one H atom in the 0-site costs about 0.2 eV to deform the 
surrounding cells if compared with absorption in the T-site. This number agrees well with the 
value previously obtained from the ab-initio simulation, 0.1 eV for a fully relaxed lattice, 
justifying the interpretation derived from these calculations. This approach shows that the 
main difference between absorption in the T-site and 0-site is related to the relaxation 
of elastic energy in neighboring cells. Only if all the cells are deformed in the same way, 
i.e. when the system does not need to invest energy again the elastic properties of the 
surrounding cells, the 0-site may become the preferred one. 



CONCLUSIONS 



For the case of high concentrations of interstitial H in BCC-iron, a large BCT distortion 
is found favorable with H absorbing preferentially in the octahedral site. This distortion is 
only possible if the hydrogen concentration is high enough to transform a significant number 
of unit cells from BCC to BCT. This transformation might form a precursor to the BCC to 
FCC phase transformation. For the case of low concentration of interstitial H, the energetic 
cost of deforming the surrounding iron matrix makes preferable absorption in the T-site. 
Hydrogen diffuses faster around regions of high concentration, where the transition state is 
located between the T-site and 0-site. In low concentration regions, diffusion takes place 
between T-site, avoiding the 0-site (a local maximum). H absorbed in 0-site or T-site create 
a completely different distribution of stresses in neighbouring cells explaining the preference 
of H for the T-site for cases or regions of low concentration of interstitials. 
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